##########################################################################################

library(data.table)
library(optparse)
library(parallel)
library(dplyr)

##########################################################################################
option_list <- list(
    make_option(c("--expression_correlation_path"), type = "character"),
    make_option(c("--tf_file"), type = "character"),
    make_option(c("--cpu"), type = "character"),
    make_option(c("--out_path"), type = "character") 
)

if(1!=1){
    
    ## 表达相关性所在的文件夹
    expression_correlation_path <- "~/20231121_singleMuti/results/tf_regulators/expression_correlation"

    ## 靶基因文件
    tf_file <- "~/20231121_singleMuti/results/tf_regulators/Positve_TF-Gene.onlyTargetGene.tsv"

    ## 线程
    cpu <- 20

    ## 输出
    out_path <- "~/20231121_singleMuti/results/tf_regulators"

}

###########################################################################################
parseobj <- OptionParser(option_list=option_list, usage = "usage: Rscript %prog [options]")
opt <- parse_args(parseobj)
print(opt)

expression_correlation_path <- opt$expression_correlation_path
tf_file <- opt$tf_file
out_path <- opt$out_path
cpu <- as.numeric(opt$cpu)

dir.create(out_path , recursive = T)

###########################################################################################

dat <- fread(tf_file)

###########################################################################################

result <- c()

for( clus in unique(dat$cluster) ){
    print(clus)

    ## 基因表达矩阵相关性
    corr_file <- paste0( expression_correlation_path , "/" , clus , "_expressionCorrelation.tsv")
    dat_corr <- fread( corr_file )
    dat_corr$FDR_expression <- p.adjust( dat_corr$pvalue , method = "fdr" )
    dat_corr$Correlation_expression <- dat_corr$correlation
    dat_corr$P_expression <- dat_corr$pvalue
    dat_corr <- dat_corr[,c("geneA" , "geneB" , "Correlation_expression" , "P_expression" , "FDR_expression")]

    ## 该细胞类型中的基因和靶基因的关系
    tmp_dat <- subset( dat , cluster == clus )

    ## TF和symbol并不和geneA、geneB对应，所以需要每行去提取表达相关系数
    result_tmp <- bind_rows(mclapply(1:nrow(tmp_dat),function(i){
        print(i)
        tmp <- tmp_dat[i,]

        ## 提取表达系数并合并
        tmp_cor <- dat_corr[which( ( dat_corr$geneA == tmp$TF & dat_corr$geneB == tmp$symbol ) | ( dat_corr$geneB == tmp$TF & dat_corr$geneA == tmp$symbol ) ),]
        tmp <- cbind( tmp , tmp_cor[,c("Correlation_expression" , "P_expression" , "FDR_expression")] )

        ## 自己和自己的表达
        if( tmp$TF == tmp$symbol ){
            tmp$Correlation_expression <- 1
        }

        return(tmp)
    },mc.cores=cpu))

    result <- bind_rows( result , result_tmp )
}

###########################################################################################
## 输出
out_file <- paste0( out_path , "/Positve_TF-Gene.onlyTargetGene.annoExpCor.tsv" ) 
write.table( result , out_file , sep = "\t" , row.names = F , quote = F )
